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A RETARDED MEAN-FIELD APPROACH FOR INTERACTING 

FIBER STRUCTURES 

R. BORSCHE*, A. KLAR*t, C. NESSLER*, A. ROTH*, AND O. TSE* 

Abstract. We consider an interacting system of one-dimensional structures modelling fibers 
with fiber-fiber interaction in a fiber lay-down process. The resulting microscopic system is inves¬ 
tigated by looking at different asymptotic limits of the corresponding stochastic model. Equations 
arising from mean-field and diffusion limits are considered. Furthermore, numerical methods for the 
stochastic system and its mean-field counterpart are discussed. A numerical comparison of solutions 
corresponding to the different scales (microscopic, mesoscopic and macroscopic) is included. 
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1. Introduction. One-dimensional structures appear in various context of in¬ 
dustrial applications. They are used, for example, in the modelling of polymers in 
suspensions, composite materials, nanostructures, fiber dynamics in turbulent flows 
and, in particular, fiber lay-down in technical processes of non-woven materials. Fur¬ 
thermore, such structures have been modelled on different levels of description involv¬ 
ing different scales. Besides microscopic models, mesoscopic kinetic or Fokker-Planck 
equations have been widely used for a statistical description of the fiber or polymer 
distributions. We refer to [1, 2] and [18, 21] for concrete examples in the industry. In 
this article, we consider non-woven materials, which are webs of long flexible fibers. 
Production processes and models corresponding to the lay-down of such fibers have 
been intensively investigated. See [4] and the above cited references. 

In the above mentioned investigations, the one-dimensional structures (fibers) un¬ 
der consideration were assumed to be mutually independent, which clearly does not 
represent reality. Therefore, the present work aims at including fiber-fiber interaction, 
thereby describing the size of each fiber and, simultaneously, the absence of intersec¬ 
tion among fibers. This is achieved by simply including the interaction of structures 
into a well investigated model described in non-woven production processes. Taking 
into account the interaction of the structures on the microscopic level leads to coupled 
systems of stochastic differential equations. Its statistical description should also take 
into account the interactions and will consequently no longer be based on the classical 
Fokker-Planck model. 

The new model makes use of a microscopic systems of retarded stochastic dif¬ 
ferential equations, and its mesoscopic description is obtained via formal mean-field 
procedures. The mean-field limit is described by a McKean-Vlasov type equation 
with a delay term. We perform an analytical investigation of the mean field limit, 
as well as a numerical comparison of microscopic, mean held and macroscopic equa¬ 
tions. The analysis of the limit is based on the work in [5, 12, 15, 16, 17, 25, 31]. For 
numerical methods for mean-held type equations we refer to [2, 26, 27, 30]. 

The paper is organized as follows: starting from a model for independent hirers, 
we present in Section 2 a new model for interacting fibers, which takes into account 
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the finite size of the fibers and prevents intersection of the fibers. The model is based 
on a system of retarded stochastic differential equations with suitable interaction po¬ 
tentials. Section 3 describes the mean-field equation and a discussion of its stationary 
solutions. The core of this section is devoted to the proof of the mean-field limit for 
the corresponding deterministic system, i.e., the rigorous derivation of the mean-field 
equation from the system of retarded deterministic equations. Section 4 contains the 
diffusion limit of the kinetic equation, while Section 5 contains a description of the 
numerical methods used for the microscopic and mean-field equations. The numerical 
results include an investigation and comparison of stationary states, as well as a con¬ 
vergence analysis of solutions to equilibrium for both the microscopic and mean-field 
equations. We finally conclude in Section 6. 

2. Interacting isotropic fiber models. We begin by reviewing a basic, mu¬ 
tually independent fiber model for the lay-down process of fibers, described by a 
stochastic dynamical system in dimensions d > 2 (cf. [4, 21, 24]). Using the state 
space A4 := x S d_1 , where S d_1 denotes the unit sphere in K d , a fiber is given by 
a path of the following stochastic differential equation: 


dx t = r t dt, 

. . / 1 \ (2.1) 

dr t = (/ - r t (g) r t ) o ^ - ——-V x V(a: t ) dt + AdW t J 

with initial condition ( xq , To) £ At. Here, V : —► ffi. is the so-called coiling potential, 
A a nonnegative diffusion coefficient and Wt the standard Brownian motion. We use 
a coordinate free formulation, where o denotes the Stratonovich integral. Note also 
that (/ — t <8> t) is the projector of onto the unit sphere S d_1 . We refer to [11] for 
related work. We equip our state space with the measure dx dv, where dx denotes the 
Lebesgue measure on R d and dv the normalized surface measure on S d_1 . Note that 
the dimensional scaling l/(d — 1) is introduced for convenience, in order to achieve a 
stationary state which does not explicitly depend on the spatial dimension. 

The fiber model above only describes the evolution of the center line of the fiber 
and does not capture the effects of the finite size of the fibers. Hence, self-intersection 
and intersection among fibers are not prevented in the model. A possibility to remedy 
this deficiency it to include hysteresis into the system, which enables the system to 
prevent self-intersection, and for multiple fibers, intersection among them. Consider¬ 
ations along this line of thought lead to the following interacting fiber model. 

Consider N £ N fibers with position and velocity (x l , r'j £ At, for each i £ 
{!,..., N}. The interacting fiber model reads 


dxl = Tf dt 

drl = {i — t\ <S> r t ‘) 


d- 1 


V x V{x\) dt 


d- 1 


i N 1 

ly 1 

N ^ t 
3=1 ' 


S7 X U ( x\ — x{) ds 


dt + A dWl 


( 2 . 2 ) 


with independent Brownian motions W t l . In comparison to the previous fiber model, 
we include a scaled, nonlocal (in time) interaction term 


V X U( a 


x J s ) ds, 


1 

t 


0 
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where U is an interaction potential, which is repulsive in our case, so as to avoid any 
contact among the fibers. Note that a ’non-retarded’ version of this system would be 
similar to a model for swarming with roosting potential (cf. [7, 8]). In addition to 
the interaction term, our new model includes a delay term, i.e., an integration with 
respect to time, to describe the fact that fibers interact with each other and with 
itself on the whole fiber length. The factor 1/N leads to the so called ’weak coupling 
scaling’ (cf. [5, 25]). Similiarly, the scaling 1/t is a normalization of the potential with 
respect to the time integration. Also note that the summation does not exclude i = j, 
which accounts for self-interaction of a fiber with itself. 

There is some freedom in the choice of the interaction potential. As a simple 
example, one can use a modifier type potential 

U(x) = U(\x\) =Cex p( - | a ,| 2 )’ for M <2R ’ 

where R is a nonnegative parameter representing the fiber radius and C > 0 is a fixed 
constant describing the strength of the interaction. Alternatively, a potential could 
be described by a smoothed version of Heaviside type potential 

U{x) = U(\x\) = CQ(2R-\x\), 


with Heaviside function 0. A smooth version of such a potential may be given by 


U{x) 


C 

1 + exp ( - k (l - (go) ) 


(2.3) 


for some regularizing parameter k > 0. 

In Figure 2.1 we compare the fiber curves for different noise amplitudes A. We 
use the coiling potential V(x) = \\x\ 2 and the interaction potential U from (2.3), with 
k = 100, C = 10 and R = 0.4. Since the computation was done for a short amount of 
time, we neglect the scaling 1 jt in front of the interaction part. We observe that the 
non-intersecting fiber curves are increasingly altered with increasing noise amplitude. 

Remark 1. 

1. For U = 0, one obtains a fully decoupled system for (x z , r l ), and each fiber is 
described by the mutually independent model given in (2.1). 

2. To consider inelastic interactions, one could include damping terms depending 
on the velocity in the interaction potential. This would lead to equations where 
a dissipative force is included in the interaction term. 

3. As in the case without interaction, it is also possible to formulate a smooth 
version of the interacting fiber system. The basic idea is to replace the Brow¬ 
nian motion on the sphere by an Ornstein-Uhlenbeck process (cf. [22]). 

Remark 2. It is also possible to include reference curves 7 into the model, which 
describes, for example, the motion of a conveyor belt. This can be done in the following 
way. We denote with rf : R + —> R d the actual fiber curves. Then (2.2) changes to 


drf t = t[ dt 

dri = (I - t\ <g> t[) 


d-1 


V x V{rf t - 7 ) dt 


d-1 


1 N 1 r * 

AM 


v x u(4 


3 l)ds 


dt + A dWl 


(2.4) 




Fig. 2.1. Comparison of the influence of the noise A for two interacting fibers. 


A change of variables £* := rf t — 7 describes the deviation of the fiber curves from the 
reference curve and thus, (2.f) may be written as 


dCt = T l t dt - dj t 

dr\ = (/ - t‘ ® r*) o ( - -±_S7 x V(fl) dt 


1 

d- 1 





£s +7t 


7 s )ds 




(2.5) 


In this case the force due to the interaction potential depends on the relative fiber point 
position as well as on the difference "ft — 7 s in the reference curve. 

In the relevant case of non-wovens on a conveyor belt moving with constant speed, 
we consider d € {2,3} and a reference curve given by "f t = -v re fe-\t. Here v re f = 
Vbeit/Vprod is the ratio between the speed of the conveyor belt, Vbeit, and the speed of the 
production process, v pro d, and e\ denotes the direction of belt movement (cf. [21]). 

In the following, we formulate a slightly more general model than (2.2). Since, 
in reality, the fiber material is transported away by the belt, interaction will not take 
place for the full history of the fibers. Thus, it is reasonable to consider a cut-off with 
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a cut-off size H > 0. Define 


t for t < H 

H for t > H ’ 


h{t) = 


He (0, oo). 


Then the interacting fiber model with cut-off is given by 


dx\ = t[ dt 

dr l t = {I -T l t ® T l t ) o ^ V x V(x l t ) dt 


1 

d- 1 


i N i 

- V — 

N h(t) 


f 

J t—h(t ) 


V x U(xl 


x{) ds 


dt + A dWl 


( 2 . 6 ) 


In the limit H — > 0, we obtain a non-retarded interacting particle model with constant 
speed, whereas, in the limit H —» oo, we obtain the interacting fiber model in (2.2). 

3. Mean-field equation. The associated mean-field equation for the distribu¬ 
tion function / = f(t,x,r ) may be formally derived from the microscopic equations 
following the procedure described, for example, in [ 6 , 17]. The equation reads 

d t f + r-\7 x f + Sf = Lf (3.1) 

with deterministic force term Sf = S v f + S u f, given by 

S V f = ~V T - (/(/-T®T) 

S u f = ~V T - ^/(J-r®r) 

where V T is the gradient on S"” 1 , and diffusion operator 

Lf = T f, (3.3) 

where A T denotes the Laplace-Beltrami operator on § d_1 . In addition, we have the 
normalization f Rd pdx = 1 , where p denotes the zeroth-moment p = f dv. 

Remark 3. A stationary solution of the mean-field equation (3.1) is given by the 
time-independent solution of 


1 


rV.R , 


d-1 

1 1 
d-1 h(t) J t _ 


/ / V x U(x-y)p(s,y)dyds , 

Jt—h(t) Jm d J 


(3.2) 


r-Vxf + Sf=^-A T f. 
Looking for a solution independent of r, we get 

p V x ^ In p + V + U * p'j = 0. 
This leads to the integral equation 

In p + V + U*p = c, 
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(3.4) 



where the constant c is fixed by the normalization f Rd pdx = 1. The integral equation 
may be written in the equivalent fixed-point form 


0 -V-U*p 


P = 


/ Rd e^' u*Pdx' 


(3.5) 


On the other hand, the stationary solution may also be characterised as the (unique) 
minimizer of the ’energy functional’ 

T(p) := / (lnp — l)pdx + / (V + U * p)pdx, 

J R d jR d 

where the first term describes the internal energy (entropy), and the second describes 
the potential energy. In contrast to the case without interaction, the stationary solu¬ 
tion has to be determined numerically. 

The rest of this section is devoted to a rigorous proof of the mean-held limit for 
the deterministic interacting fiber model 


dxj 

dt 

dr[ 

dt 




1 T, O - 


1 




d- 1 

1 1 

d - 1 N h(t) J t -h(t) 


(3-6) 


1 


\7 X U(x\ - x 3 s ) ds , 


towards the Vlasov type equation 

dtf + r■ V x f + Sf = 0, 
where S = S v + S u is as given in (3.2). 


(3.7) 


3.1. Mean-field limit of a retarded system. We consider a system of N £ N 
interacting particles with state Z\ £ R m , m £ N at time t £ R + 


dZ( . 1 \ ^ 1 

-rr = a(Z t ) + hJI) 


h (t) Jt- 


h(t) 


B(Zi,Zi)ds, Zq = z i £ 


(3.8) 


H £ (0,oo), or h(t) = t, 


for each i £ { 1,..., N}, where as above 

, ft for t < H 
m= \H for OH' 

and a £ Lip(R m ), B £ Lip fc (R m x R m ;R m ) be globally Lipscliitz, satisfying 

sup \B(zuz) - B(z 2 ,z)\ < Lip(U) 1 A \zi - z 2 \, 

z£R m 

sup \B(z,z 1 ) - B(z,z 2 ) | < Lip (B) 1 A\z\ — z 2 \, 
ze R™ 

where we use the notation x/\y = min{a;, y}. As in (2.2) and (2.6), (3.8) is a retarded 
system of ordinary differential equations. 
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Denote by V r (R m ), r > 0, the set of Borel probability measures p such that 

[ \z\ r p(dz) < oo. 

J R m 

For every probability measure p £ Vi(R m ), we set 

JCp(z) := p(B(z, •)) = [ B(z,z)p(dz), z G R m . (3.9) 

J R m 

Taking the empirical measure p^ = $(' ~ Z() e V\(W n ) as p in (3.9) gives 

N 

^s(Z l t) = / B{Z' t ,x)tf{dx) = -Y,B(Zi,Zi). 

Jr™ j=i 

Therefore, we may consider a mean-field equation given by 

dZ 1 /'* 

-jri z ) = F (t,Z t (z),{fj,} t ):=a(Z t (z)) +jj-r K.p s (Z t (z)) ds, (3.10) 

at h(t) J t _ h ( t) 

with initial condition Zq(z) = z, law (Zq) = pg. Here fit = Z t (-)#po denotes the push 
forward of the measure po € i.e., 

f tp(z) ix t (dz) = ( ip(Z t (z)) po(dz) \/ipeC b (M. m ), 

Jr™ Jr™ 

under the flow Z £ C(M + x R m ;R m ) generated by the mean-field equation, and {p}t 
denotes the family of measures {p s , s £ [0, t]} up to time t > 0. 

We begin with an existence and uniqueness result for solutions of (3.10). Our 
results generalizes those given in [16], which follows the general scheme introduced 
in [12] (cf. [5, 25]). The proof relies on a slight modification of the proof in [16, 
Proposition 4.1]. For the sake of completeness, we include the proof in Appendix A. 

Proposition 3.1. Let the assumptions on a and B above be satisfied. Then the 
mean-field equation (3.10) admits a unique global solution Z £ C(R+ x R m ;R m ). 

Remark 4. The family of measures {pt = Z t (-)#po} c Ri(K m ) generated by the 
flow Z £ C(R_|_ x R m ;R m ) provides a weak solution to the Vlasov equation 

d t pt + div(F(t,z, {p}t)pt) = 0, (3.11) 

More precisely, for all h £ C“(R m ), the functions pt(h) are differentiable, 

^M=p t (F(t,.,{fx} t ).Vh) 

for all t > 0, and pt(h) —> po(h) for t —> 0+. 

Next, we show that solutions to the mean-field equation (3.10) depend continu¬ 
ously on the initial probability measures po G Ri(R m ). To do so, we need to measure 
the difference of two probability measures p, v £ P 1 (R m ). A convenient way is to use 
the Monge-Kantorovich-Rubinstein distance W\ on Ri(]R m ) defined in [20] (cf. [32]), 


Wi(p,v) = inf 

7TGII 


1 A \x — y\ ir(dxdy), 
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where II(/i, v) is the set of Borel probability measures on R m x K m such that 



{<t>{x) + ip(y))n{dx dy) 


(j>(z) dz + 


il>{z) dz 


for all G C&(R m ). This distance is also called the Wasserstein distance. 

Proposition 3.2. Let p J 0 G V\ (K m ), and Zi G C(K+ x R m ;R m ) be the corre¬ 
sponding solution to the mean-field equation (3.10) with initial conditions 

Z 1 0 (z) = z, law(Zg) = p 3 0 , for j = 1,2 


then for every T > 0, the following stability estimate holds 
sup < c(l + T)e cT 

t€[0,T] 


for some constant c = c(F, H) > 0, where p J t := Z\ 

Proof. Let ttq G II(/xJ,/Zo), and define 

V(t) ■■= [[ 1 A \Z$(zi) - Z((z 2 )\TT 0 (dz 1 dz 2 ) 

J «/R m xl m 

Following the arguments in the proof of Proposition 3.1, we derive the estimate 

T>(t) < 2?(0) + Lip(F) / V(s) ds + Lip(£?) [ —-j— [ V(a)dads. 

Jo Jo 'H s ) Js-h(s) 

We only show the estimate for the case H = oo, i.e. h(t) = t. The general case 
H > 0 may be shown analogously. Similar to Proposition 3.1, by simple but tedious 
computations, we have 

V(t) < D(0) + Lip(F) f g(t,s)V(s)ds, 

Jo 

with g(t, s) = 1 + ln(t) — ln(s). Recursively, we obtain 


1 - Lip(F) 


fc (n + l )T h 


k\ 


sup V(t) < VLip(P) 
te[o,T] 


/(* + l )T l 


i-o 


l\ 


Therefore, passing to the limit k —> oo yields 


sup V(t) < c(l+T)e cT 2?(0), 

i€[0,T] 


with c = max{l,Lip(T 1 )}. Since 

m = 


1 A \zx - z 2 \n t (dz 1 dz 2 ), 


V(0). 


where ir t is the push forward measure of tt 0 under the map (Z), Zf) 1 and 


7T 0 g n (pl,pl) 


n t G n 



for any t > 0, we have 


sup W\ (nl, p() < c(l + T)e cT V(0) 
te[o,T] 

Finally, optimizing in ttq yields the required assertion. □ 

Remark 5. For any N G N, the family of empirical measures {p^,t > 0} C 
'Pi(R m ) defined by fi^ = Z t (-)#fi(( G Vi(R m ), where 

1 N 

~ ^ 3 = 1.JV 

j =i 

is a weak solution to the Vlasov equation (3.11) by construction. Therefore, if we 
know that limjv->-cx> Pq) —> 0 for some /io G Vi(R m ), then the stability result 

given in Proposition 3.2 provides the convergence 

lim Wiint,^) -> 0 for any t G [0,T], 

A—>oo 

where pt = Z t (-)#ii o, with Z G C([0,T] x R 777- ;!! 171 ). 

3.2. Application to the retarded fiber equations. We now use the results 
above to show the mean-field limit of (3.6) towards (3.7). For this reason, we denote 
Z\ = (; x\,t\) G R d x R d , and write a = (a\ ,a 2 ), B = (B 1 ,B 2 ) with 

ai (^ ; ) = tI a 2 (Zl) = -J—(l- T t® T () o V x V(xl), 

Bi(Zi, Zi) = 0, B 2 (Zl Zi) = V x U(xi - xi). 

Obviously, we are unable to directly apply the results developed above, since a and 
B do not satisfy the assumptions above. Consider B 2 for the moment. Tedious but 
simple computations yield 

\B 2 (z 1 ,z) - B 2 (z 2 ,z)\ < Lip(V x £/)((l + |v-!| 2 )|a?i - x 2 \ + (|ti| + |t 2 |)|ti - r 2 |) 
\B 2 (z,zi) - B 2 (z,z 2 )\ < Lip(V x t/)(l + |r| 2 )|xi - x 2 \ 

Hence, B is not globally Lipscliitz. Fortunately, if we only consider B on the subset 
M C K d x R d , then Lip(R) = X U). Consequently, B is globally Lipschitz on 

AT Clearly, the same conclusion holds for a. 

In order to ensure that Z\ G M. for all t > 0, we observe that 

^l T t| 2 = 0 for all t > 0, i = 1,... ,N. 

Therefore, if we start with r* G S d_1 , then also t) G S d_1 for all other times t > 0, and 
we may apply our results obtained above to initial measures /io G "Pi(R d x R d ) with 
supp(/io) C AI. Indeed, since Z t (x,r) G A! for any ( x,t ) G M, we may consider the 
flow on AI and the push forward /i t = Z t (-)#p .o G 'Pi(AI), as soon as supp(/io) C AT 
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4. Large diffusion scaling. In this section we formally investigate situations 
with large values for the noise amplitude on a diffusive time scale, i.e., we change 
L = eL and t = et (cf. [3, 19]). Replacing the scaled operators in (3.1) and omitting 
the tilde lead to the scaled equation 


edtf + r- V x f + Sf = -Lf. 

e 


We use a Hilbert expansion of the form / = f 0 + efi 
simply get /o = fo{x) = p(x). To order e, one obtains 


+ ... for (4.1). To order 1, we 


T-V x f 0 + Sfo=—&rfl, 


which, due to the Fredholm alternative, gives 


h = - 


A 2 (d- 1) 


r-foVx ln/o + V+ — 


h{t) Jt-h 


{U * h){s,-)ds . 


Integrating (4.1) with respect to dv gives 


fdv + V a 


rf dv = 0. 


Considering terms up to order e, we obtain 


d t fo + Vz • / t/i dv = 0. 

J s^- 1 


Therefore, inserting f\ and computing the integral over the tensor product yields 


{'" l, + v + W)L m (u * f ' ){a '' )ds ) ■ 

This equation is similar to an equation derived in [4, 19, 24] for the case without an 
interaction potential. The stationary equation reads 

pV x { lnp + V + U * p\ = 0, / pdx = 1, 

v ' Jr* 

which leads again to the integral equation (1.1) as for the mean-field case. 

5. Numerical Method and Results. In this section we investigate the quali¬ 
tative behavior of solutions corresponding to the interacting fiber model numerically. 
More specifically, we consider the case of isotropic fibers in three spatial dimensions. 

5.1. Numerical methods. We describe numerical solvers for the microscopic, 
mean-field and macroscopic equations, respectively. 

5.1.1. Microscopic model. To solve equations (2.2) or (2.6) numerically we 
use the Euler Maruyama method. Using Ito integration (2.2) is written as 


dx\ = rf dt 


dri = -\{l-Tl ® r*) (\7 x V{x\) + V x U{x\ - x{) ds^j dt 


N ^ t 

3 = 1 



+ A 2 t l dt + A(I - t\ <g> Tl)d,W l t . 


The Ito form of (2.6) is obtained in an analogous way. We consider an equidistant 
time grid given by 0 = t 0 < ... < t n with step size At. We denote x l n := x\ for the 
position of fiber i at time t n = nAt. The time integration in (2.2) is approximated by 


1 

t 


V x U(x l t 


x{) ds 


E ~ < )At = l E - <)■ (5-i) 

; 1 n ; , 


For (2.6) with H > 0, we set up a buffer to store the previous positions x\. The buffer 
size is determined by nf, u f = h(t)/At. The time integration is then approximated by 

i rt nbuf i 

— / V x U(x\ - xi) ds ~ E - <_ fc )At. (5.2) 

' l W Jt—h(t) n \ t n) 

Hence, the Euler-Maruyama iteration is given by 


X n+1 X r 


T n+1 ~ T « 


At-rl 
At ■ 


■^Y,(Y,T W toU{x i n ~xl)At)) -A 2 4 


N 


i 


V x V(xi) 


N ■ 


t, 


3=1 fc=l 

+ VAt-A(l-r^4)K 


where R^ £ M 3 , for i = 1 is a vector containing three normally distributed 

pseudorandom values. For (2.6), one has to replace the sum over all previous time 
steps (5.1) by the sum over the time step buffer (5.2). Typical simulations of the 
microscopic system include N ~ 10 6 fibers, whose positions are used to obtain a 
histogram, which will be compared to the solution of the mean field equation described 
below. For such large values of TV, the evaluation of the interaction term is quite costly. 
Several measures were taken to ensure reasonable computing time: 

1. We did not sum up every time step in (5.1) and considered only every n-th 
time step. In that way, the results are not altered in a significant way, if h is 
not too large. 

2. Furthermore, the Fortran code carrying out the microscopic simulations was 
parallelized using OpenMP. If a reasonable scaling of the computing time with 
the number of processors is to be achieved, the parts of the code running in 
parallel have to be as mutually independent as possible, in order to reduce 
the overhead in communication between parallel threads. 

3. Therefore, we divided N into smaller groups of fibers, which then interact 
only within each group, thereby enabling a parallel computation with very 
little overhead. The histogram for the spatial density is then produced from 
the position data of all the groups. 

5.1.2. Mean field equation. The numerical methods used here are an advance¬ 
ment of the schemes used in [9, 27, 28]. We apply a second order Strang splitting [10] 
to (3.1) to obtain subproblems on spatial and velocity domain. We split the equation 


dtf + T- Z x f + Sf = Lf, 
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Fig. 5.1. The left figure depicts two connected triangles of the geodesic grid, while the right 
figure shows the full spherical grid. 


into the subequations 


dtf {1) = ~\(Sf W -Lf^), (5.3) 

d t f {2) = -T-V x f^\ (5.4) 

dtf {3) = ~\(SfW-LfW). (5.5) 

Equation (5.4) has to be solved only on the spatial domain R 3 , while (5.3) and (5.5) 

are to be solved on the velocity domain § 2 . Note however, that all the computation 
steps on one grid have to be carried out for every grid point on the other grid. 

Firstly, we discuss the solution of (5.3) and (5.5) with a finite volume discretization 
of the velocity space based on a geodesic grid (cf. [27, 29]). The geodesic grid consists 
of spherical triangles, which are represented with vertices and normals in R 3 , neighbor 
relations and correct distance and surface measures. In this way the geometry is 
included implicitly into the method. We denote by Ti the i-th cell of the grid with 
cell midpoint t* . The cell midpoint is chosen to be the intersecting point of the great 
circle arcs passing perpendicularly through the midpoints of the cell edges. | ■ | denotes 
a length or surface measure, T i3 denotes the edge between cells T* and X), 7y ; denotes 
the edge midpoint, and e(r,j) denotes the outer normal vector of cell i at the edge 
midpoint r- i; . The distance from cell midpoint r, to Tj is given by hi ,, which is divided 
into the parts h\ and / 12 . Due to the construction of the grid, the distance between 
cell midpoints does not vary that much, and hi, hi ~ h^/2 , so all the cells have nearly 
the same size and shape. The reader is referred to Figure 5.1 for a visualization of 
the grid structure and notations. 

As for any finite volume scheme, the solution is obtained via cell averages 


/" = 7^7 J t f{t„> x,t) dr, 

where dr is the canonical surface measure on the sphere. We define 

F(t,x,r) := \{I - t ® t)(\/ x V + 7-^— [ \7 x U*pds ), 

^ V fl \T) Jt-h(t) ' 
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and integrate (5.3, 5.5) over the cells [t n ,t n+ 1 ] x T,; to obtain 


rn+1 _ rn _ ^ \ 

Ji Ji - | T .| 


j£N(i) 


F(s, x , r) • e(r)/ dr ds 


+ 


rtn+i 


V r / • e(r) dr ds 


Applying the midpoint rule on the cell edge and in time, one obtains the iteration 

1 


f? +1 f? 

At 


mi 


E 

j€N(i) 


I Tij\F(x, Tij ) ■ e{Tij)f(t n+ 1 , x, Tij ) 


A 2 


The overall order of the method depends on the approximation of f(t n+ i,x,Tjj) and 
the normal flux V r /(t n+ i, x, Tij) ■ e(rjj). In this case, we choose 


f{t n+ i,X,Tij) 


hi + ~y{F{x, nj) ■ e(Tij)) h 2 - -£{F{x, t tj ) ■ e( Tij )) 

u J i ' i J ?' 


Vr./'(/„ . i :X.Tij) ■ e(jij) 


f? ~ f? 

hi. i 


The approximation for f(t n+ ±,x,Tij) is similar to the Lax-Wendroff numerical flux 
function, except that we evaluate F(x,t) ■ e(r) at the edge midpoint Tjj, and not at 
the cell midpoints Tj and Tj, since the grid structure only provides normal vectors 
at the cell interfaces. The value for V T /(i n+ i, x, Tjj) ■ e(rjj) is obtained by a finite 
difference approximation on the connecting circle arc of the cell midpoints r,; and Tj , 
see [27]. Although the method is not a second order method, the numerical results 
are close to those of a second order method, see the discussion below and Figure 5.2 
for the convergence rates of the splitting scheme. 

Equation (5.4) is solved using a semi-Lagrangian method [23, 30] on an equidistant 
grid Xjjk G M 3 ,i,j, k = 1, ...,n x with grid size Ax. The characteristic curves q(t) of 
(5.4) starting at grid point Xjjk at time t n+ 1 are given by 


7 (t) = Xjjk + t-T. 


The derivative of / with respect to time on a characteristic curve is given by 
7 (t),T) = d t f(y(t), T,t)+T ■ V x f(j(t), T , t)= 0 
and it follows, that / is constant along the characteristic curve: 

f{t n +i,Xjjk,T) = /(t„,7(-At),r). 

Since only the values of / at the grid points at time t n are known, we have to 
interpolate f(t n , y(—At), r) from f(t n ,Xjjk,T). The order of the method depends 
on the order of the interpolation procedure, since the characteristic curves can be 
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computed analytically in this case. An order higher than one will produce unphys¬ 
ical oscillations at discontinuities of the numerical solution, which has to be pre¬ 
vented. We use a Bezier interpolation [9]. We use the notations Xijk = (x\j k )i— 1 , 2 , 3 , 
y(-At) = (7i(- Ai ))z=i,2,3 e [x} jk ,xl jk + Ax] x [xf jk ,x? jk + Ax] x [x? jk ,xf jk + Ax] 
and define £ = (£i)i= 1 , 2,3 as follows 


6 


li (~At) - x l ijk 
Ax 


l = 1,2,3. 


Then the interpolated values on each cell are given by an interpolating polynomial of 
three space variables of the following form: 


3 

f(t n ,l(-At),T)= ^2 S Aj3 (5 1 )B At; 3(^2)S 1/j3 (^3)x A/ily , 

A, / 1 , i/=0 


where H; i3 (£) = ( 3 )£*(1 — ^ € {0,1,2,3} are the cubic Bernstein polynomials, 

and i> A m „ are the control values for the interpolation. This interpolant never leaves 
the convex hull of the control values. These values are chosen appropriately, such 
that the required interpolation order is preserved on smooth sections of the solution, 
and oscillations are prevented elsewhere. This can be done by computing an interpo¬ 
lating Newton polynomial of appropriate order, performing a change of basis to the 
Bernstein polynomials and then shifting control points back into the convex hull of 
the neighboring grid function values. See [9] for more information on this method 

We note that semi-Lagrangian methods are not conservative in general. For lin¬ 
ear advection, we get a conservative method using the Bezier interpolation procedure 
without slope limiting. In [26], this is shown for a third order interpolating polyno¬ 
mial, which is identical to the Newton polynomial. Since we reproduce this polynomial 
in the Bernstein polynomial basis, the same computation can be done for the Bezier 
interpolation used here. However, the slope limiting procedure destroys mass conser¬ 
vation and we have to apply conservation techniques as described in [14, 23]. 

5.2. Stationary equation. The stationary equation (3.4) is solved via an iter¬ 
ation scheme. We use the fixed-point form (3.5) and the iterative scheme 

e -{V+U*p n ) 

Pn+1 = fe~( v + u *P")dx' 

As starting point p 0 we use the solution of the stationary equation without interaction, 
namely po = Ce~ l . The iteration is well-defined, i.e., p n is integrable and has integral 
one. In addition, p n is strictly positive for all n € N. 

For the implementation we have to approximate the convolution U*p or S7 X U* p. 
This is done using the midpoint rule on each grid cell. One obtains 

(U*p)(x ijk )= U(x ijk - y)p(y)dy « ^ U(x ijk - x AjLtl/ )p(x A ^)(Ax) 3 , (5.6) 


where x^k are the grid points and Ax is the (constant) grid size in each direction, 
as before. For the computation of the convolution, one has to compute the distance 
matrix (| x^k — x AAiy |), which consumes a lot of computing time and memory. One 
advantage over the solution of the microscopic system is, that the grid points are fixed 
in contrast to the fiber positions in the microscopic context. Therefore, the distance 
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matrix may be precomputed. Furthermore, one can precompute for each grid point 
a list of neighbor points, which have a relevant effect on the computation for a given 
set of parameters and the purely repulsive potential (2.3). Only these points provide 
a relevant contribution to the interaction force. 

The current implementation only considers interaction forces, which are higher 
than 0.1 percent of the maximal occurring force as relevant. A considerable amount 
of memory can be saved by applying this weak restriction. The time integration in 
(3.2) for h(t) = t does not increase the effort in the mean-field context, because it can 
be realized by just summing up the interaction forces from different time steps. For 
the case H G (0, oo), a time step buffer as in the microscopic context is used. 

5.3. Numerical results. We show numerical results for the microscopic model 
(2.6) and the mean-field equation (3.1) in three spatial dimensions. In particular, 
we investigate the convergence of the microscopic and mean-field solutions to the 
corresponding stationary state. The stationary states obtained from microscopic and 
mean-field computations are compared to the solution of (3.5). For a more detailed 
investigation of the convergence to equilibrium for different values of the delay ff, we 
investigate the time decay of the distance to equilibrium for microscopic and mean- 
field equations and compare it to the one obtained for the case without interaction. 
Moreover, the sensitivity of the stationary solutions with respect to the parameters 
of the interaction model is investigated. 

For all computations we consider the three dimensional case for (2.6). Unless 
otherwise stated, we choose the interaction potential given in (2.3) with k = 10, 
C = 10 and R = 1.4. The coiling potential V is chosen as V(x) = \x\ 2 /2. 

We are especially interested in the comparison of the numerical results for the 
mean-held equation and the microscopic stochastic model. For this comparison, we 
consider different noise coefficients A and a box-shaped initial condition for the mean- 
held equation given by 


for x G [—1, l] 3 , r 3 > 0 
else 

where r 3 denotes the third component of the vector r and C is a normalizing con¬ 
stant. The corresponding initial condition for the microscopic system can be obtained 
by randomly placing particles in the cube [—1, l] 3 with initial velocities having the 
third component larger than zero. We start the investigation by studying order of 
convergence of the numerical method for the mean-held equation described above. 

5.3.1. Order of convergence of the numerical method for the mean-field 
equations without interaction. In this section we give a numerical investigation 
of the method for the mean held equation for the case without interaction. Although 
the method is not a second order method, owing to the geodesic grid, the numerical 
results show that the numerical order of convergence is of order two. Let n x be the 
number of spatial grid points in each direction and n*, the number of cells on the 
sphere, h x and hk are the respective grid sizes. Note, that hk is the average grid size 
on the geodesic grid. The spatial step size was chosen to match hk- The following 
grid sizes were simulated: 


fo(x,r) =C 
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convergence 



Fig. 5.2. Order of convergence for the mean-field equation without interaction. The green and 
red lines correspond to first and second order error, respectively. 


n x 

11 

21 

41 

n k 

20 

80 

320 

hx 

0.76 

0.38 

0.19 

hk 

0.730 

0.353 

0.175 

err L 2 

0.004604 

0.000777 

0.000175 


From this data, we obtain the convergence plot shown in Figure 5.2. 

5.3.2. Stationary solution. In Figure 5.3, we consider the stationary distribu¬ 
tions obtained from the microscopic and mean-field equations for different values of 
the delay H. We show radial plots of the radially symmetric density function. The 
case without interaction is compared to the case with interaction, with delay given 
by H £ {0, 0.1, 0.5, oo}. The grid resolutions simulated for the mean-field equations 
were n x = 40 points in each spatial direction with a grid size of h x = 0.205 for the 
cases without interaction and the cases H £ {0, oo}. To slightly reduce the computing 
time, we used a smaller spatial resolution of n x = 35 and h x = 0.235 in the cases with 
H £ {0.1,0.5}. In all cases, the sphere was discretized by 320 geodesic triangles, which 
led to an average cell midpoint distance of 0.175. These grid sizes where matched by 
the histograms generated from the microscopic fiber positions and velocities. 

To achieve smooth data, the total number of simulated fibers has to be chosen 
sufficiently large. In the case without interaction, the microscopic density is generated 
from 4.8 • 10 6 simulated fiber positions. For the cases with interaction, we computed 
800 realizations of a system of 600 interacting fibres and generated the histogram 
from the fibre positions of all the realizations, which means that the histogram is 
based on 4.8 • 10 5 microscopic fiber positions. As mentioned earlier, the interaction 
computation in the microscopic context is quite costly. So, although the computing 
time for 4.8 • 10 6 fibres without interaction took only a couple of minutes, depending 
on how many processing cores we have at our disposal, it took up to several days 
of computing time for the case H = 0.5, which is the most expensive task here. 
Computing times for the mean-field solver lie in the range of 3 to 24 hours on the 
given hardware, which was an Intel XEON E5 2670 with 8 cores at 3.3GHz. 

As expected, the stationary distribution does not change with different values of 
H. However the stationary distribution for the case with interaction is significantly 
wider than the distribution for the case without interaction. We note that the solu- 
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mean field, stationary 


SDE, stationary 




Fig. 5.3. In blue: stationary solution without interaction; in green: stationary solution for 
retarded interaction term without cut-off, in red: stationary solution for retarded interaction term 
with cut-off H = 0.1. 


tions obtained by the microscopic and by the mean-field equation are in very good 
agreement for the cases considered here. 

5.3.3. Convergence to equilibrium. In this subsection we invstigate the con¬ 
vergence to equilibrium in more detail. In Figures 5.4 and 5.5, we investigate the 
behaviour (in time) of the L 2 -norm of the distance to equilibrium for the spatial 
density p. Again, the case without interaction is compared to the case with inter¬ 
action with cut-off given by H £ {0,0.1, 0.5, oo}. In Figure 5.4, the case A = 1.0 is 
considered, while Figure 5.5 shows the case A = 0.5. 


dist. to equilibrium mean field, A=1.0 


dist. to equilibrium SDE, A=1.0 




Fig. 5.4. In blue: decay without interaction; in green: decay for retarded interaction term 
without cut-off, in red: decay for retarded interaction term with different H and A = 1.0. 


The figures provide a strong indication that all cases with interaction converge to 
the same stationary distribution. We also observe the exponential decay to equilibrium 
for the case without interaction (cf. [13] for a theoretical exposition). Furthermore, 
the graph in green shows the case with a strongly retarded potential. The decay is 
no longer strictly monotone and it takes more time to achieve the stationary state, in 
comparison to the cases without significant delay. 

Figure 5.6 displays a comparison between the decay to equilibrium for different 
values of A. As can be observed, the decay becomes slower and, in particular, less 
and less regular for smaller values of A. 
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dist. to equilibrium mean field, A=0.5 


dist. to equilibrium SDE, A=0.5 




Fig. 5.5. In blue: decay without interaction; in green: decay for retarded interaction term 
without cut-off, in red: decay for retarded interaction term with different H and A = 0.5. 



Fig. 5.6. Comparison of decay to equilibrium for different values of A with H = oo. 

5.3.4. Sensitivity of stationary solutions with respect to parameters in 
the interaction potential. Finally, we investigate the behavior of the stationary 
solution for different parameters in the interaction model. In particular, we expect 
a wider profile for larger radii, as well as for stronger interaction, i.e., C large. To 
numerically analyze the behavior for the different parameters, we consider the sta¬ 
tionary equation in fixed-point form (3.5). We are interested in the change of the 
stationary solution for varying parameters R and C respectively. 

In Figure 5.7, the stationary solution is shown. First, we consider the case of 
varying radius R. We use the potential U from (2.3) and keep C and k fixed. As 
we can see in Figure 5.7, the solution p is becoming less concentrated for increasing 
radius R. Second we vary C keeping R and k fixed obtaining similiar results as before. 

6. Conclusion and outlook. In this work we extended existing models for 
describing the lay-down process of fibers during the production of non-wovens to a 
model, which captures the interaction between the fibers by using a system of retarded 
stochastic differential equations and their corresponding mean-field approximations. 
A theoretical investigation of the mean-field limit for the deterministic retarded sys¬ 
tem is included. The solutions of the (retarded) microscopic and mean-field approx¬ 
imations are numerically compared with each other. In particular, the stationary 
solutions and the decay to equilibrium is numerically investigated. Different types 
of interaction potentials are presented and the influence of different parameters in 
the interaction potential is investigated. Furthermore, a large diffusion scaling for 
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Fig. 5.7. Stationary solutions p for different values of the radius R and different values of the 
factor C in the interaction potential. 


the mean-field equations has been considered. Another interesting theoretical issue 
which has not been treated in the present work is the analytical investigation of the 
convergence of the retarded mean-field equation towards equilibrium. 

Appendix A. Proof of Proposition 3.1. Here we give a proof of Proposition 
3.1, compare [16]. We begin by writing the differential equation in its integral form 

Z t (z)=z + [ a(Z s (z)) ds + [ — F- [ JC/j, a (Z s (z))dads. (1.1) 
JO Jo 'H s ) Js-h(s) 

Let T > 0 be arbitrary but fixed and set Bt '■= C([0, T\ x R' m ;K' m ). We now denote 
T to be the operator on the right hand side, and show the existence of a unique fixed 
point for this operator. 

To see that this operator is well-defined, let Z £ Bt- The first to terms on the right 
hand side are obviously well-defined. Let us now consider the third term. By definition 
of the push forward measure, n = Z#/i 0 £ C([0,T];w-Pi(K m )). Furthermore, since 

\ICnt(zi) - K.nt(z 2 )\ < [ \B(zi,z) - B(z 2 ,z)\ m(dz) < Lip(B)|zi - z 2 \, 

A" 


for every t £ [0,T], we infer that /C/r £ C([0,T]; Lip b (R m )) ^ Bt- A simple reformu¬ 
lation of the term under the integral leads to 


1 


h (s) J s - 


h(s) 


^/^<r(^s(•2’)) d(J — / K-'fta(s,h(s),cr) (^s (-^)) d(J 

Jo 

n B{Z s {z), Z a ( sMs)>a) (z)) Vo{dz) da , 

i m 


with a(s, h(s), a) = (1 — a)(s — h(s)) + as. From the assumptions, we conclude that 
this term is continuous in s £ [0, T\ and Lipschitz in z £ K m . Consequently, the third 
term is well-defined. Altogether, we have that T: Bt —► Bt as required. 

As a first estimate, we obtain 


\(TZ-TZ) t (z)\<Up(F) / \Z s (z)-Z s (z)\ds 


+ Lip(B) f f ( | Z a {z) 
Jo ”'( s ) Js-h(s) Jr™ 


Z a (z) | no(dz) da ds, 
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from the Lipschitz continuity of the mappings on the right hand side. Denoting 
£(t) = sup \{TZ — TZ) t {z)\ and £(t) := sup \(Z — Z) t (z)\, 

zgR m zSR™ 

we obtain for the cases H G (0, oo), t < H and h(t) = t. the following estimate: 

£(t) < Lip(.F') [ £(s) ds + Lip(B) [ - [ £(a)dads 

Jo Jo s Jo 

= Lip(F) f £(s) ds + Lip(i3) f (ln(t) — ln(s)) £(s) ds 
Jo Jo 

< Lip(F) f (1 + ln(t) — ln(s)) £(s) ds, 

Jo 

where we used integration by parts in the equality. With this estimate at hand, we 
construct a solution by the Picard iteration procedure, and show that the generated 
sequence converges. We, therefore, define the sequence (Z^)fc>o C £>t recursively by 
Z (k+i) = 7 ~z^) for k > 0, and denote 

£ k (t) = sup |(Z (fc+1) - Z (k) ) t {z) |, g(t,s) = 1 + ln(<) - ln(s). 

z6R™ 

Applying the above estimate iteratively leads to 

£ k (t)<Lip(F) f g(t,ti)£ k -i(t 1 )dt 1 
Jo 

< Lip(f 1 ) 2 f g(t,t i) f g{ti,t 2 ) £k- 2 {t- 2 ) dt 2 dti 
Jo Jo 


< Lip(^) A 


rt rtk-i 

/ g(t, ti) ■ s / g(t k -i,t k ) dt k ■ sdti 
Jo Jo 


sup £o{tk). 

0 <t k <t 


Elementary computations of the term in the bracket gives 

£k(t) < Lip (F) k ^ k + ^ t sup £ 0 (t k ), 
K- o <t k <t 


for any k > 0. Since the series 

^ {k + 1 )t^ , v C £ 

2 c k -— —- = (1 + ct)e 4 < +oo, 

k> o 

summing up the terms £k{t) in k > 0 yields 

V ]£k(t) < sup £ 0 (t k )(l + Up(F)t)e hip{F)t < +oo, 

k> 0 0<t k <t 

and hence £ k (t) —> 0 pointwise in t > 0. Consequently, the Picard sequence (Z ( ' k ' , ) k > 0 
converges uniformly in Bt to the solution of the integral equation (1.1). Since T was 
chosen arbitrarily, the solution may be extended to Z £ C(M+ x R m ;R m ). 
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As for uniqueness, we take two solutions Z and Z of the integral equation, and 
denote its difference by 


£(t) = sup |(Z - Z) t (z )|. 

zg K m 


From the estimates above, we have 

£(t) < Lip(F) f g(t, s) £(s) ds. 

Jo 

Let T > 0 be arbitrary but fixed. Following the arguments from above, we obtain 


1 - Lip(F) 


fc ( k + l)t k 

k\ 


sup £(t) < 0. 

o <t<T 


Passing to the limit k — > oo yields £ = 0 on [0, T\, and hence Z = Z on [0,T] x R m . 
Since T was chosen arbitrarily, this uniqueness result extends to R+ x R m . 

We now consider the case H G (0, oo). As above, we establish a unique solution 
Y € Br- Setting Yh(z) = y as initial value, the integral form for t > H reads 

z t{y) = y+ [ a(Z s (y))ds+ f -[ Ky a (Z s (y)) da ds. 

JH JH H s ) Js-h(s) 

With the same notations as above, we obtain the following estimate 


£(t) < Lip(F) / £ (s) ds + Lip(i3) 
J H 


1 


£(cr) da ds. 


Ih H s ) Js-h(s) 

The second term on the right hand side may be expressed in the form 

[ [ £(a)dads=— f £(s)ds— [ £(s)ds 

JH his) Js-h(s) H Jt-H Jo 

- f H s (£( s ) -£(s-H))ds 

t r t r H 

= — / £(s)ds—l £(s)ds 

H Jt-H Jo 


1 [ T 1 r*- 11 

- jjJ H S£( s )ds+ — (s + H)£(s) 

t /■* 

~ H Jo £ ^ ds ’ 

where we used integration by parts in the first equality, and the fact that 
rt pt — H pt—H pt 

/ {H — s)£(s') ds < 0, / (s — t)£(s) ds < 0, / £(s)ds< / £{s)ds , 

Jh Jo Jo Jo 

in the last inequality. Consequently, there is a constant c = c(F, H) > 0, such that 


ds 


£(t)<c(l+t) / £(s)ds. 
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Constructing a sequence (Z^)k >o £ Bt via the Picard iteration for some arbitrary 
but fixed T > H, with zf 1 = Y t for t £ [0, H], we obtain as above 


£k(t) < (1 +t)c k 


t k t k+1 \ 
k\ + (k + 1 )!) 


sup £ 0 (t k ). 

o <t k <t 


( 1 . 2 ) 


Notice that £k = 0 on [0 ,H], Replicating the arguments above yields existence and 
uniqueness of a solution Z £ Bt- Moreover, we may extend the solution to M+ x R m . 
To conclude, we observe that 


f [ B(Z( k \z),zil MsM (z))^dz)da 

J 0 J R 771 

— > B(Z s (z),Z a{sMs) ^(z)) ^ 0 (dz)da, 

Jo Jr™ 


for k —> oo by the Lebesgue dominated convergence, and so 

i r 


h(s) J s _ 


h(s) 


JC/j, a (Z s (z)) da 


is continuous for all z € R m . Consequently, t Z t (z), t > 0, is continuously differ¬ 
entiable for all z £ R m , and satisfies the differential form (3.10). □ 
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